Compiled under R version 3.6.2 (2019-12-12)
This document details all analyses performed in R for the study: Legendre, L. J., Rubilar-Rogers, D., Musser, G. M., Davis, S. N., Otero, R., Vargas, A. O., and Clarke, J. A. 2020. A giant soft-shelled egg from the Late Cretaceous of Antarctica. Nature.
For more information regarding the study, datasets, and analyses, please refer to the Supplementary Information of this paper. If you have any additional questions, feel free to email me at lucasjlegendre@gmail.com.
library(AICcmodavg)
library(ape)
library(caper)
## Loading required package: MASS
## Loading required package: mvtnorm
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
##
## select
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(evobiR)
## Registered S3 method overwritten by 'geiger':
## method from
## unique.multiPhylo ape
##
## Attaching package: 'evobiR'
## The following object is masked from 'package:AICcmodavg':
##
## AICc
library(ggplot2)
library(MPSEM)
library(nlme)
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
library(nortest)
library(phylopath)
library(phytools)
## Loading required package: maps
library(RColorBrewer)
This part requires the use of dataset 1 (“Dataset1-lepidosaurs.txt”; see also Supplementary Table 1 in Legendre et al., 2020), and lepidosaur phylogenetic tree (“Lepidosaurtree.trees.nex”).
tree<-read.nexus("Lepidosaurtree.trees.nex")
data2<-read.table("Dataset1-lepidosaurs.txt", header=T); data2V<-subset(data2, !is.na(V)&!is.na(SVL))
pruned.tree<-drop.tip(tree, setdiff(tree$tip.label, data2V$Species))
data<-read.table("Dataset1-lepidosaurs.txt", header=T, row.names = "Species"); dataV<-subset(data, !is.na(V)&!is.na(SVL))
alpha <- seq(0, 1, 0.1)
fit <- list()
form <- log(V)~log(SVL)
for (i in seq_along(alpha)) {
cor <- corMartins(alpha[i], phy = pruned.tree, fixed = T)
fit[[i]] <- gls(form, correlation = cor, data = data, na.action=na.exclude, method = "ML")
}
sapply(fit, logLik)
## [1] Inf -319.2756 -323.6809 -324.3353 -324.6167 -324.8171 -324.9815
## [8] -325.1226 -325.2464 -325.3564 -325.4550
BM<-gls(log(V)~log(SVL), data=dataV, correlation=corBrownian(phy=pruned.tree), method="ML")
OU<-gls(log(V)~log(SVL), data=dataV, correlation=corMartins(1, phy=pruned.tree, fixed=T), method="ML")
Lambda<-gls(log(V)~log(SVL), data=dataV, correlation=corPagel(1, phy=pruned.tree), method="ML")
EB<-gls(log(V)~log(SVL), data=dataV, correlation=corBlomberg(0.1, phy=pruned.tree), method="ML")
OLS<-gls(log(V)~log(SVL), data=dataV, method="ML")
Cand.models = list()
Cand.models[[1]] = BM
Cand.models[[2]] = OU
Cand.models[[3]] = Lambda
Cand.models[[4]] = EB
Cand.models[[5]] = OLS
Modnames = paste(c("BM", "OU", "Lambda", "EB", "OLS"), sep = " ")
aictab(cand.set = Cand.models, modnames = Modnames, sort = T)
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## Lambda 4 485.22 0.00 1 1 -238.52
## EB 4 550.19 64.97 0 1 -271.01
## BM 3 571.55 86.33 0 1 -282.72
## OU 3 657.01 171.79 0 1 -325.45
## OLS 3 658.92 173.70 0 1 -326.41
summary(Lambda)
## Generalized least squares fit by maximum likelihood
## Model: log(V) ~ log(SVL)
## Data: dataV
## AIC BIC logLik
## 485.0497 498.9889 -238.5248
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.8802218
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.7990941 0.7306413 -1.093689 0.2752
## log(SVL) 1.6708400 0.0941323 17.749917 0.0000
##
## Correlation:
## (Intr)
## log(SVL) -0.63
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.11666418 -0.80674408 -0.01263292 0.65520753 1.74914328
##
## Residual standard error: 1.168844
## Degrees of freedom: 241 total; 239 residual
plot(Lambda, resid(., type="n")~fitted(.), main="Normalized Residuals v Fitted Values",abline=c(0,0))
res<-resid(Lambda, type="n")
par(mar=c(5.1,4.1,4.1,2.1))
qqnorm(res)
qqline(res)
lillie.test(residuals(Lambda))
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuals(Lambda)
## D = 0.07408, p-value = 0.002707
lillie.test(chol(solve(vcv(pruned.tree)))%*%residuals(Lambda))
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: chol(solve(vcv(pruned.tree))) %*% residuals(Lambda)
## D = 0.15828, p-value = 4.371e-16
which(res>4); which(res<(-2.5))
## Cryptoblepharus_poecilopleurus
## 168
## Chirindia_ewerbecki Micrurus_mipartitus
## 89 222
newdata<-dataV[-c(89,168,222),]
newdata2<-data2V[-c(89,168,222),]
newtree<-drop.tip(pruned.tree, setdiff(pruned.tree$tip.label, newdata2$Species))
Newlambda<-gls(log(V)~log(SVL), data=newdata, correlation=corPagel(1, phy=newtree), method="ML")
summary(Newlambda)
## Generalized least squares fit by maximum likelihood
## Model: log(V) ~ log(SVL)
## Data: newdata
## AIC BIC logLik
## 428.773 442.6621 -210.3865
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.9221827
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.8303776 0.7046654 -1.1784 0.2398
## log(SVL) 1.6870012 0.0865058 19.5016 0.0000
##
## Correlation:
## (Intr)
## log(SVL) -0.6
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.24045710 -0.84905671 -0.03920481 0.62244401 1.55314354
##
## Residual standard error: 1.1406
## Degrees of freedom: 238 total; 236 residual
plot(Newlambda, resid(., type="n")~fitted(.), main="Normalized Residuals v Fitted Values",abline=c(0,0))
res<-resid(Newlambda, type="n")
par(mar=c(5.1,4.1,4.1,2.1))
qqnorm(res)
qqline(res)
lillie.test(residuals(Newlambda))
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuals(Newlambda)
## D = 0.074716, p-value = 0.002609
lillie.test(chol(solve(vcv(newtree)))%*%residuals(Newlambda))
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: chol(solve(vcv(newtree))) %*% residuals(Newlambda)
## D = 0.13226, p-value = 7.181e-11
#newlambda:
datacomp<-comparative.data(phy=newtree, data=newdata2, names.col="Species")
pgls<-pgls(log(V)~log(SVL), data=datacomp, lambda="ML")
summary(pgls)
##
## Call:
## pgls(formula = log(V) ~ log(SVL), data = datacomp, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.107910 -0.039203 -0.005097 0.034535 0.147102
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.718
## lower bound : 0.000, p = 1.8499e-05
## upper bound : 1.000, p = 4.5223e-05
## 95.0% CI : (0.413, 0.908)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.41450 0.91690 -0.4521 0.6528
## log(SVL) 1.63215 0.15212 10.7293 6.661e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06154 on 64 degrees of freedom
## Multiple R-squared: 0.6427, Adjusted R-squared: 0.6371
## F-statistic: 115.1 on 1 and 64 DF, p-value: 6.121e-16
#lambda:
datacompl<-comparative.data(phy=pruned.tree, data=data2V, names.col="Species")
pgls<-pgls(log(V)~log(SVL), data=datacompl, lambda="ML")
summary(pgls)
##
## Call:
## pgls(formula = log(V) ~ log(SVL), data = datacompl, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.140589 -0.048249 -0.004924 0.031043 0.155446
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.715
## lower bound : 0.000, p = 7.7636e-05
## upper bound : 1.000, p = 0.00015352
## 95.0% CI : (0.381, 0.918)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.19647 0.96534 0.2035 0.8393
## log(SVL) 1.53068 0.16093 9.5115 5.418e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06683 on 66 degrees of freedom
## Multiple R-squared: 0.5782, Adjusted R-squared: 0.5718
## F-statistic: 90.47 on 1 and 66 DF, p-value: 5.417e-14
Not much difference between the two models, assumptions of normality are not met in either of them – we use lambda instead of newlambda for the plot.
ggplot(dataV, aes(log(SVL), log(V), color=Suborder)) +
geom_point(size=4) +
xlab("ln snout-vent length (mm)") +
ylab("ln egg volume (mm3)") +
geom_abline(intercept=-0.7990941, slope=1.6708400, colour="lightblue", size=1.3) +
theme(panel.background = element_rect(fill="black")) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_colour_brewer("Suborder", palette="Accent")
exp((-0.7990941/1.6708400)+(1/1.6708400)*log(5471405.789))
## [1] 6684.303
data[1,7]<-6684.303
dataV2<-subset(data, !is.na(V))
ggplot(dataV2, aes(log(SVL), log(V), color=Suborder)) +
geom_point(size=4) +
xlab("ln snout-vent length (mm)") +
ylab("ln egg volume (mm3)") +
geom_abline(intercept=-0.7990941, slope=1.6708400, colour="lightblue") +
theme(panel.background = element_rect(fill="black")) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_colour_brewer("Suborder", palette="Accent")
Eggshell thickness ~ body mass is not detailed here, since the code is virtually identical – only variable names and regression coefficients for the plot and estimation are different.
LambdaBM<-gls(log(V)~log(BM), data=dataV, correlation=corPagel(1, phy=pruned.tree), method="ML")
summary(LambdaBM)
## Generalized least squares fit by maximum likelihood
## Model: log(V) ~ log(BM)
## Data: dataV
## AIC BIC logLik
## 466.1175 480.0567 -229.0587
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.7332279
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 5.050556 0.4341909 11.63211 0
## log(BM) 0.619536 0.0300611 20.60926 0
##
## Correlation:
## (Intr)
## log(BM) -0.254
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.1068229 -0.3986971 0.2864704 0.8395741 2.8852589
##
## Residual standard error: 0.9261633
## Degrees of freedom: 241 total; 239 residual
exp(-(5.050557/0.619536)+(1/0.619536)*log(5471405.789))
## [1] 21657214
data[1,8]<-21657214
dataV2<-subset(data, !is.na(V))
ggplot(dataV2, aes(log(BM), log(V), color=Suborder)) +
geom_point(size=4) +
xlab("ln body mass (g)") +
ylab("ln egg volume (mm3)") +
geom_abline(intercept=5.050556, slope=0.619536, colour="lightblue") +
theme(panel.background = element_rect(fill="black")) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_colour_brewer("Suborder", palette="Set1")
dataVT2<-subset(data2, !is.na(Thickness)&!is.na(V)); dataVT2<-dataVT2[-1,]
pruned.treeVT<-drop.tip(tree, setdiff(tree$tip.label, dataVT2$Species))
dataVT<-subset(data, !is.na(Thickness)&!is.na(V)); dataVT<-dataVT[-1,]
BMVT<-gls(log(Thickness)~log(V), data=dataVT, correlation=corBrownian(phy=pruned.treeVT), method="ML")
OUVT<-gls(log(Thickness)~log(V), data=dataVT, correlation=corMartins(1, phy=pruned.treeVT), method="ML")
LambdaVT<-gls(log(Thickness)~log(V), data=dataVT, correlation=corPagel(1, phy=pruned.treeVT), method="ML")
EBVT<-gls(log(Thickness)~log(V), data=dataVT, correlation=corBlomberg(0.1, phy=pruned.treeVT, fixed=T), method="ML")
OLSVT<-gls(log(Thickness)~log(V), data=dataVT, method="ML")
Cand.models = list()
Cand.models[[1]] = BMVT
Cand.models[[2]] = OUVT
Cand.models[[3]] = LambdaVT
Cand.models[[4]] = EBVT
Cand.models[[5]] = OLSVT
Modnames = paste(c("BMVT", "OUVT", "LambdaVT", "EBVT", "OLSVT"), sep = " ")
aictab(cand.set = Cand.models, modnames = Modnames, sort = T)
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## EBVT 3 134.92 0.00 0.73 0.73 -64.27
## BMVT 3 137.75 2.83 0.18 0.91 -65.69
## LambdaVT 4 139.22 4.30 0.09 1.00 -65.29
## OLSVT 3 148.59 13.67 0.00 1.00 -71.11
## OUVT 4 150.85 15.93 0.00 1.00 -71.11
summary(EBVT)
## Generalized least squares fit by maximum likelihood
## Model: log(Thickness) ~ log(V)
## Data: dataVT
## AIC BIC logLik
## 134.5443 141.2028 -64.27216
##
## Correlation Structure: corBlomberg
## Formula: ~1
## Parameter estimate(s):
## g
## 0.1
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 2.0900026 0.4535976 4.607614 0
## log(V) 0.3270484 0.0574324 5.694490 0
##
## Correlation:
## (Intr)
## log(V) -0.972
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.9317947 -0.7662902 -0.2506952 0.4323433 3.0468614
##
## Residual standard error: 0.6624366
## Degrees of freedom: 68 total; 66 residual
lillie.test(residuals(EBVT))
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuals(EBVT)
## D = 0.073029, p-value = 0.4927
lillie.test(chol(solve(vcv(pruned.treeVT)))%*%residuals(EBVT))
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: chol(solve(vcv(pruned.treeVT))) %*% residuals(EBVT)
## D = 0.12937, p-value = 0.006569
plot(EBVT, resid(., type="n")~fitted(.), main="Normalized Residuals v Fitted Values",abline=c(0,0))
res<-resid(EBVT, type="n")
par(mar=c(5.1,4.1,4.1,2.1))
qqnorm(res)
qqline(res)
which(res<(-2)); which(res>2.5)
## Zootoca_vivipara_viviparous
## 45
## Chamaeleo_senegalensis
## 19
newdataVT<-dataVT[-c(19,45),]
newdataVT2<-dataVT2[-c(19,45),]
newtreeVT<-drop.tip(pruned.treeVT, c("Zootoca_vivipara_viviparous","Chamaeleo_senegalensis"))
NewEBVT<-gls(log(Thickness)~log(V), data=newdataVT, correlation=corBlomberg(0.1, phy=newtreeVT, fixed=T), method="ML")
summary(NewEBVT)
## Generalized least squares fit by maximum likelihood
## Model: log(Thickness) ~ log(V)
## Data: newdataVT
## AIC BIC logLik
## 114.2644 120.8334 -54.13222
##
## Correlation Structure: corBlomberg
## Formula: ~1
## Parameter estimate(s):
## g
## 0.1
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 2.1246786 0.4052994 5.242244 0
## log(V) 0.3209485 0.0511337 6.276657 0
##
## Correlation:
## (Intr)
## log(V) -0.972
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.7567958 -0.8478322 -0.2694998 0.4934623 2.2910465
##
## Residual standard error: 0.5849121
## Degrees of freedom: 66 total; 64 residual
plot(NewEBVT, resid(., type="n")~fitted(.), main="Normalized Residuals v Fitted Values",abline=c(0,0))
res<-resid(NewEBVT, type="n")
par(mar=c(5.1,4.1,4.1,2.1))
qqnorm(res)
qqline(res)
lillie.test(residuals(NewEBVT))
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuals(NewEBVT)
## D = 0.085226, p-value = 0.2753
lillie.test(chol(solve(vcv(newtreeVT)))%*%residuals(NewEBVT))
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: chol(solve(vcv(newtreeVT))) %*% residuals(NewEBVT)
## D = 0.12101, p-value = 0.01773
NewEBVT$coefficients
## (Intercept) log(V)
## 2.1246786 0.3209485
ggplot(newdataVT, aes(log(V), log(Thickness), color=Suborder)) +
geom_point(size=4) +
xlab("ln egg volume (mm3)") +
ylab("ln eggshell thickness (µm)") +
geom_abline(intercept=2.1246786, slope=0.3209485, colour="lightblue") +
theme(panel.background = element_rect(fill="black")) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_colour_brewer("Suborder", palette="Accent")
ppadata<-dataVT
ppadata[,c(4,6:8)]<-log(ppadata[,c(4,6:8)])
M<-define_model_set(null=c(),
one=c(Thickness~SVL),
two=c(V~BM),
three=c(Thickness~V),
four=c(Thickness~V+SVL),
five=c(Thickness~V, V~BM),
six=c(Thickness~SVL, V~BM),
seven=c(Thickness~SVL+BM),
eight=c(V~BM, Thickness~BM),
nine=c(Thickness~V+BM),
ten=c(Thickness~SVL+BM+V),
eleven=c(Thickness~V+BM, V~BM),
twelve=c(Thickness~SVL+BM, V~BM),
thirteen=c(Thickness~BM),
.common=c(BM~SVL, V~SVL))
plot_model_set(M)
## Warning: The curvature argument has been deprecated in favour of strength
result<-phylo_path(M, data=ppadata, tree=pruned.treeVT, model='lambda'); result
## A phylogenetic path analysis, on the variables:
## Continuous: Thickness SVL V BM
## Binary:
##
## Evaluated for these models: null one two three four five six seven eight nine ten eleven twelve thirteen
##
## Containing 31 phylogenetic regressions, of which 10 unique
s<-summary(result); s
## model k q C p CICc delta_CICc l w
## eleven eleven 1 9 0.212 8.99e-01 21.3 0.000 1.00e+00 5.13e-01
## five five 2 8 2.982 5.61e-01 21.4 0.108 9.48e-01 4.86e-01
## eight eight 2 8 20.246 4.46e-04 38.7 17.372 1.69e-04 8.67e-05
## nine nine 2 8 21.093 3.04e-04 39.5 18.218 1.11e-04 5.68e-05
## three three 3 7 23.863 5.53e-04 39.7 18.415 1.00e-04 5.15e-05
## four four 2 8 21.871 2.13e-04 40.3 18.996 7.50e-05 3.85e-05
## twelve twelve 1 9 20.009 4.52e-05 41.1 19.797 5.03e-05 2.58e-05
## six six 2 8 23.369 1.07e-04 41.8 20.494 3.55e-05 1.82e-05
## ten ten 1 9 20.881 2.92e-05 42.0 20.669 3.25e-05 1.67e-05
## two two 3 7 35.727 3.12e-06 51.6 30.278 2.66e-07 1.37e-07
## thirteen thirteen 3 7 41.127 2.73e-07 57.0 35.679 1.79e-08 9.18e-09
## seven seven 2 8 40.890 2.83e-08 59.3 38.015 5.56e-09 2.85e-09
## one one 3 7 45.578 3.59e-08 61.4 40.129 1.93e-09 9.92e-10
## null null 4 6 57.935 1.18e-09 71.3 49.997 1.39e-11 7.14e-12
plot(s)
## Warning: Position guide is perpendicular to the intended axis. Did you mean to
## specify a different guide `position`?
## Warning: guide_axis(): Discarding guide on merge. Do you have more than one
## guide with the same position?
## Warning: guide_axis(): Discarding guide on merge. Do you have more than one
## guide with the same position?
## Warning: guide_axis(): Discarding guide on merge. Do you have more than one
## guide with the same position?
## Warning: guide_axis(): Discarding guide on merge. Do you have more than one
## guide with the same position?
averagemodel<-average(result)
## Registered S3 methods overwritten by 'MuMIn':
## method from
## nobs.pgls caper
## nobs.phylolm phylolm
## logLik.phylolm phylolm
plot(averagemodel)
## Warning: The curvature argument has been deprecated in favour of strength
fossiltreePEM<-read.nexus("Lepidosaurtree.trees.nex")
fossildataPEM<-read.table("Dataset1-lepidosaurs.txt",header=TRUE,stringsAsFactor=FALSE,)
fossildataV<-subset(fossildataPEM, !is.na(V))
fossildataV[,c(3:5,7:9)]<-log(fossildataV[,c(3:5,7:9)])
treePEM<-drop.tip(fossiltreePEM, setdiff(fossiltreePEM$tip.label, fossildataV$Species))
treedrop <- drop.tip(treePEM,fossildataV[is.na(fossildataV[,"SVL"]),"Species"])
grloc <- getGraphLocations(treePEM,fossildataV[is.na(fossildataV[,"SVL"]),"Species"])
sporder <- match(attr(grloc$x,"vlabel")[grloc$x$vertex$species],fossildataV[,"Species"])
PEMfs <- list()
PEMfs[["V"]] <- PEM.fitSimple(y=fossildataV[sporder,"SVL"],
x=fossildataV[sporder,"V"],w=grloc$x,d="distance",sp="species",
lower=0,upper=1)
PEMfs[["none"]] <- PEM.fitSimple(y=fossildataV[sporder,"SVL"],
x=NULL,w=grloc$x,d="distance",sp="species",lower=0,upper=1)
for(m in c("V","none")) print(PEMfs[[m]]$optim$par)
## [1] 0.4981507
## [1] 0.1197493
rm(m)
PEMAIC <- list()
PEMAIC[["V"]] <- lmforwardsequentialAICc(y=fossildataV[sporder,"SVL"],
x=fossildataV[sporder,"V",drop=FALSE],object=PEMfs[["V"]])
PEMAIC[["none"]] <- lmforwardsequentialAICc(y=fossildataV[sporder,"SVL"],object=PEMfs[["none"]])
for(m in c("V","none"))
cat(m,summary(PEMAIC[[m]])$adj,PEMAIC[[m]]$AICc,"\n")
## V 0.9945104 -248.9384
## none 0.9913577 -124.1351
rm(m)
summary(PEMAIC[["V"]])
##
## Call:
## lm(formula = as.formula(paste(p1, p2, sep = "")), data = df1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.178636 -0.035423 0.004214 0.039172 0.153809
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.071635 0.066712 31.053 < 2e-16 ***
## V 0.393493 0.009214 42.705 < 2e-16 ***
## V_1 -9.100852 0.108193 -84.116 < 2e-16 ***
## V_2 -4.704946 0.126827 -37.097 < 2e-16 ***
## V_30 1.372985 0.083670 16.410 < 2e-16 ***
## V_6 0.986086 0.110409 8.931 3.32e-15 ***
## V_10 0.858355 0.118694 7.232 3.56e-11 ***
## V_8 -0.984514 0.083734 -11.758 < 2e-16 ***
## V_63 -0.760769 0.084320 -9.022 1.98e-15 ***
## V_7 0.862697 0.083710 10.306 < 2e-16 ***
## V_31 0.874262 0.084017 10.406 < 2e-16 ***
## V_117 0.819127 0.084435 9.701 < 2e-16 ***
## V_27 -1.240447 0.106374 -11.661 < 2e-16 ***
## V_77 -0.733867 0.083703 -8.768 8.30e-15 ***
## V_229 -0.819382 0.086023 -9.525 < 2e-16 ***
## V_156 0.507389 0.085388 5.942 2.39e-08 ***
## V_124 -0.408268 0.089180 -4.578 1.08e-05 ***
## V_5 0.388503 0.087559 4.437 1.91e-05 ***
## V_180 -0.572372 0.083671 -6.841 2.72e-10 ***
## V_14 -0.450315 0.085040 -5.295 4.87e-07 ***
## V_225 -0.546724 0.083745 -6.528 1.33e-09 ***
## V_48 -0.617529 0.083862 -7.364 1.77e-11 ***
## V_11 0.632745 0.084051 7.528 7.39e-12 ***
## V_108 0.675610 0.084757 7.971 6.76e-13 ***
## V_76 0.648948 0.084352 7.693 3.05e-12 ***
## V_44 -0.384803 0.085978 -4.476 1.64e-05 ***
## V_17 0.485665 0.083976 5.783 5.10e-08 ***
## V_80 0.535035 0.083675 6.394 2.61e-09 ***
## V_205 0.551436 0.083744 6.585 1.00e-09 ***
## V_23 -0.401879 0.084055 -4.781 4.61e-06 ***
## V_59 -0.320047 0.085352 -3.750 0.000265 ***
## V_66 0.544710 0.084217 6.468 1.81e-09 ***
## V_13 -0.499332 0.083867 -5.954 2.26e-08 ***
## V_29 -0.594629 0.085593 -6.947 1.57e-10 ***
## V_170 0.403913 0.083791 4.820 3.90e-06 ***
## V_109 -0.541708 0.084637 -6.400 2.54e-09 ***
## V_98 -0.558359 0.085053 -6.565 1.11e-09 ***
## V_179 -0.390290 0.083857 -4.654 7.86e-06 ***
## V_18 -0.529632 0.085127 -6.222 6.15e-09 ***
## V_206 -0.439153 0.083739 -5.244 6.13e-07 ***
## V_54 -0.453173 0.083883 -5.402 3.00e-07 ***
## V_125 -0.458915 0.084018 -5.462 2.28e-07 ***
## V_142 0.402447 0.083669 4.810 4.08e-06 ***
## V_235 0.429425 0.083804 5.124 1.05e-06 ***
## V_155 0.370097 0.083716 4.421 2.04e-05 ***
## V_34 0.351834 0.083814 4.198 4.94e-05 ***
## V_194 0.273009 0.085159 3.206 0.001692 **
## V_183 -0.376601 0.083670 -4.501 1.48e-05 ***
## V_169 -0.344993 0.083766 -4.119 6.71e-05 ***
## V_112 0.278725 0.084640 3.293 0.001274 **
## V_39 -0.293362 0.084304 -3.480 0.000682 ***
## V_55 -0.391209 0.083709 -4.673 7.26e-06 ***
## V_69 0.380311 0.083680 4.545 1.24e-05 ***
## V_129 -0.395590 0.083754 -4.723 5.89e-06 ***
## V_173 -0.381593 0.083692 -4.559 1.16e-05 ***
## V_152 -0.131500 0.088366 -1.488 0.139119
## V_56 0.147425 0.087294 1.689 0.093629 .
## V_62 0.424582 0.084164 5.045 1.48e-06 ***
## V_104 0.367887 0.083752 4.393 2.29e-05 ***
## V_72 -0.370855 0.083789 -4.426 2.00e-05 ***
## V_88 -0.321197 0.083667 -3.839 0.000191 ***
## V_178 0.381534 0.084033 4.540 1.26e-05 ***
## V_133 -0.435337 0.084932 -5.126 1.04e-06 ***
## V_58 0.361402 0.083871 4.309 3.19e-05 ***
## V_68 0.363490 0.084099 4.322 3.03e-05 ***
## V_172 -0.275813 0.083706 -3.295 0.001266 **
## V_176 0.264259 0.083676 3.158 0.001972 **
## V_151 0.262899 0.083679 3.142 0.002077 **
## V_101 -0.279064 0.083671 -3.335 0.001109 **
## V_130 0.354518 0.084439 4.199 4.93e-05 ***
## V_20 0.450383 0.087594 5.142 9.68e-07 ***
## V_16 0.364896 0.084960 4.295 3.38e-05 ***
## V_4 0.448939 0.088214 5.089 1.22e-06 ***
## V_134 -0.322932 0.084261 -3.833 0.000196 ***
## V_45 -0.265568 0.083671 -3.174 0.001874 **
## V_120 0.273162 0.083703 3.263 0.001404 **
## V_227 -0.273278 0.083759 -3.263 0.001408 **
## V_33 0.347907 0.085649 4.062 8.33e-05 ***
## V_148 0.267991 0.083752 3.200 0.001725 **
## V_118 -0.283059 0.083948 -3.372 0.000982 ***
## V_138 0.267968 0.083770 3.199 0.001731 **
## V_121 -0.253752 0.083678 -3.032 0.002925 **
## V_42 0.299744 0.084538 3.546 0.000544 ***
## V_97 -0.272391 0.083932 -3.245 0.001489 **
## V_116 -0.260448 0.083813 -3.107 0.002314 **
## V_73 -0.209929 0.083934 -2.501 0.013613 *
## V_85 0.228299 0.083680 2.728 0.007242 **
## V_19 0.264901 0.084013 3.153 0.002003 **
## V_28 -0.307979 0.086077 -3.578 0.000486 ***
## V_26 0.224630 0.083682 2.684 0.008207 **
## V_21 -0.292578 0.085859 -3.408 0.000871 ***
## V_61 -0.242060 0.083830 -2.888 0.004544 **
## V_110 -0.242602 0.083878 -2.892 0.004479 **
## V_177 0.245034 0.084113 2.913 0.004208 **
## V_193 -0.215603 0.083688 -2.576 0.011095 *
## V_99 -0.244967 0.085178 -2.876 0.004704 **
## V_95 -0.219343 0.084007 -2.611 0.010080 *
## V_187 0.201047 0.083676 2.403 0.017677 *
## V_115 -0.194921 0.083706 -2.329 0.021409 *
## V_40 -0.204034 0.083748 -2.436 0.016183 *
## V_35 0.231018 0.086910 2.658 0.008837 **
## V_38 0.223027 0.085936 2.595 0.010530 *
## V_150 -0.203677 0.084038 -2.424 0.016732 *
## V_222 -0.198031 0.083780 -2.364 0.019564 *
## V_168 -0.188724 0.083980 -2.247 0.026296 *
## V_141 0.183297 0.083777 2.188 0.030449 *
## V_12 0.180644 0.083706 2.158 0.032745 *
## V_91 -0.174506 0.084069 -2.076 0.039873 *
## V_201 -0.168706 0.083691 -2.016 0.045863 *
## V_67 -0.162550 0.083667 -1.943 0.054182 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08367 on 131 degrees of freedom
## Multiple R-squared: 0.997, Adjusted R-squared: 0.9945
## F-statistic: 399.9 on 109 and 131 DF, p-value: < 2.2e-16
The model with V as a co-predictor has the best fit.
tf<-list(V=identity, none=identity)
m <- "V" ; atr <- tf[[m]](fossildataV[is.na(fossildataV[,"SVL"]),m,drop=FALSE])
resultsPEM<-predict(object=PEMfs[[m]],targets=grloc,lmobject=PEMAIC[[m]],newdata=atr,interval="confidence")
# Predicted SVL, with upper and lower limits of the confidence interval for the prediction:
exp(resultsPEM$values); exp(resultsPEM$upper); exp(resultsPEM$lower)
## [,1]
## Antarcticoolithus_bradyi 6665.867
## [,1]
## Antarcticoolithus_bradyi 7888.606
## [,1]
## Antarcticoolithus_bradyi 5632.654
The code to predict missing values for body mass is identical to the one for SVL; simply replace ‘SVL’ with ‘BM’.
data[1,7]<-6665.867
dataV2<-subset(data, !is.na(V))
ggplot(dataV2, aes(log(SVL), log(V), color=Suborder)) +
geom_point(size=4) +
xlab("ln snout-vent length (mm)") +
ylab("ln egg volume (mm3)") +
geom_abline(intercept=-0.7990941, slope=1.6708400, colour="lightblue") +
theme(panel.background = element_rect(fill="black")) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_colour_brewer("Suborder", palette="Accent")
data[1,8]<-456096.7
dataV2<-subset(data, !is.na(V))
ggplot(dataV2, aes(log(BM), log(V), color=Suborder)) +
geom_point(size=4) +
xlab("ln body mass (g)") +
ylab("ln egg volume (mm3)") +
geom_abline(intercept=5.050556, slope=0.619536, colour="lightblue") +
theme(panel.background = element_rect(fill="black")) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_colour_brewer("Suborder", palette="Set1")
RT<-data$Thickness/data$L; data<-cbind(data, RT)
dataA<-subset(data, !is.na(RT)&!is.na(OV))
dataA$RT<-log(dataA$RT)
anovatree<-drop.tip(tree, setdiff(tree$tip.label, rownames(dataA)))
datANOVA<-as.data.frame(cbind(dataA$RT, dataA$OV), header=T)
row.names(datANOVA)<-row.names(dataA); colnames(datANOVA)<-c("RT", "OV")
datANOVA<-ReorderData(anovatree, datANOVA, taxa.names="row names")
phylANOVA(anovatree, datANOVA$OV, datANOVA$RT, p.adj="BH")
## Warning: no labels for x. Assuming order of tree$tip.label.
##
## Warning: no labels for y. Assuming order of tree$tip.label.
## ANOVA table: Phylogenetic ANOVA
##
## Response: y
## Sum Sq Mean Sq F value Pr(>F)
## x 2.884385 2.884385 6.074072 0.006
## Residual 31.816186 0.474868
##
## P-value based on simulation.
## ---------
##
## Pairwise posthoc test using method = "BH"
##
## Pairwise t-values:
## 1 2
## 1 0.000000 2.464563
## 2 -2.464563 0.000000
##
## Pairwise corrected P-values:
## 1 2
## 1 1.000 0.006
## 2 0.006 1.000
## ---------
ggplot(dataA, aes(x=OV, y=RT, fill=OV)) +
geom_boxplot() +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_fill_hue()
log(data$RT)[1] # Value for Antarcticoolithus, intermediate between the groups
## [1] 0.9120565
p2<-ggplot(dataA, aes(x=OV, y=RT, fill=OV)) + geom_boxplot()
ggplot_build(p2)$data
## [[1]]
## fill ymin lower middle upper ymax outliers
## 1 #F8766D -0.1251631 0.9808293 1.35720692 1.7429693 2.784423 3.442019
## 2 #00BFC4 -0.3462762 -0.3308263 -0.02868824 0.8270896 2.503459
## notchupper notchlower x flipped_aes PANEL group ymin_final ymax_final xmin
## 1 1.5065672 1.2078466 1 FALSE 1 1 -0.1251631 3.442019 0.625
## 2 0.8860653 -0.9434418 2 FALSE 1 2 -0.3462762 2.503459 1.625
## xmax xid newx new_width weight colour size alpha shape linetype
## 1 1.375 1 1 0.75 1 grey20 0.5 NA 19 solid
## 2 2.375 2 2 0.75 1 grey20 0.5 NA 19 solid
datANOVA2<-as.data.frame(cbind(dataA$RT, dataA$Strategy), header=T)
row.names(datANOVA2)<-row.names(dataA); colnames(datANOVA2)<-c("RT", "Strategy")
datANOVA2<-ReorderData(anovatree, datANOVA2, taxa.names="row names")
phylANOVA(anovatree, datANOVA2$Strategy, datANOVA2$RT, p.adj="BH")
## Warning: no labels for x. Assuming order of tree$tip.label.
##
## Warning: no labels for y. Assuming order of tree$tip.label.
## ANOVA table: Phylogenetic ANOVA
##
## Response: y
## Sum Sq Mean Sq F value Pr(>F)
## x 7.86154 3.930770 9.666176 0.007
## Residual 26.83903 0.406652
##
## P-value based on simulation.
## ---------
##
## Pairwise posthoc test using method = "BH"
##
## Pairwise t-values:
## 1 2 3
## 1 0.000000 -1.534310 3.960211
## 2 1.534310 0.000000 4.382854
## 3 -3.960211 -4.382854 0.000000
##
## Pairwise corrected P-values:
## 1 2 3
## 1 1.0000 0.3380 0.0015
## 2 0.3380 1.0000 0.0015
## 3 0.0015 0.0015 1.0000
## ---------
ggplot(dataA, aes(x=Strategy, y=RT, fill=Strategy)) +
geom_boxplot() +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_fill_hue()
p2<-ggplot(dataA, aes(x=Strategy, y=RT, fill=Strategy)) + geom_boxplot()
ggplot_build(p2)$data
## [[1]]
## fill ymin lower middle upper ymax outliers
## 1 #F8766D -0.1251631 0.9794977 1.3139737 1.74296931 2.7844232
## 2 #00BA38 0.3137979 1.3231783 1.6757606 2.06667501 2.5034588 3.442019
## 3 #619CFF -0.3462762 -0.3359763 -0.3256763 -0.02868824 0.2682998
## notchupper notchlower x flipped_aes PANEL group ymin_final ymax_final xmin
## 1 1.47662916 1.1513182 1 FALSE 1 1 -0.1251631 2.7844232 0.625
## 2 2.02995340 1.3215677 2 FALSE 1 2 0.3137979 3.4420194 1.625
## 3 -0.04536403 -0.6059886 3 FALSE 1 3 -0.3462762 0.2682998 2.625
## xmax xid newx new_width weight colour size alpha shape linetype
## 1 1.375 1 1 0.75 1 grey20 0.5 NA 19 solid
## 2 2.375 2 2 0.75 1 grey20 0.5 NA 19 solid
## 3 3.375 3 3 0.75 1 grey20 0.5 NA 19 solid
This part requires the use of dataset 2 (“Dataset2-amniotes.txt”; see also Supplementary Table 3 in Legendre et al., 2020), and amniote phylogenetic tree (“Amniotetree.nex”).
dataS<-read.table("Dataset2-amniotes.txt", header=T)
treeS<-read.nexus("Amniotetree.nex")
treeS<-drop.tip(treeS, setdiff(treeS$tip.label, dataS$Taxon))
dataS<-read.table("Dataset2-amniotes.txt", header=T, row.names="Taxon")
Ws<-diag(vcv.phylo(treeS))
dataS<-read.table("Dataset2-amniotes.txt", header=T)
datahard<-dataS[-c(50,86:88,93:95,115,118,120:145,147,149),]
treehard<-drop.tip(treeS, setdiff(treeS$tip.label, datahard$Taxon))
dataS<-read.table("Dataset2-amniotes.txt", header=T, row.names="Taxon")
datahard<-dataS[-c(50,86:88,93:95,115,118,120:145,147,149),]
Wh<-diag(vcv.phylo(treehard))
alpha <- seq(0, 1, 0.1)
fit <- list()
form <- log(Eggshell_thickness)~log(Egg_mass)
for (i in seq_along(alpha)) {
cor <- corMartins(alpha[i], phy = treehard, fixed = T)
fit[[i]] <- gls(form, correlation = cor, data = datahard, na.action=na.exclude, weights=varFixed(~Wh), method = "ML")
}
sapply(fit, logLik)
## [1] Inf -34.53566 -35.81116 -36.82220 -37.35043 -37.58546 -37.68186
## [8] -37.72003 -37.73494 -37.74073 -37.74298
g <- seq(0.1, 1, 0.1)
fit <- list()
form <- log(Eggshell_thickness)~log(Egg_mass)
for (i in seq_along(g)) {
cor <- corBlomberg(g[i], phy = treehard, fixed = T)
fit[[i]] <- gls(form, correlation = cor, data = datahard, na.action=na.exclude, weights=varFixed(~Wh), method = "ML")
}
sapply(fit, logLik)
## [1] -47.38547 -58.17006 -63.69614 -67.03638 -69.29137 -70.92493 -72.16651
## [8] -73.14384 -73.93423 -74.58743
BM<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datahard, correlation=corBrownian(phy=treehard), weights=varFixed(~Wh), method="ML")
OU<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datahard, correlation=corMartins(0.1, phy=treehard, fixed=T), weights=varFixed(~Wh), method="ML")
Lambda<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datahard, correlation=corPagel(1, phy=treehard), weights=varFixed(~Wh), method="ML")
EB<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datahard, correlation=corBlomberg(0.1, phy=treehard, fixed=T), weights=varFixed(~Wh), method="ML")
OLS<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datahard, method="ML")
Cand.models = list()
Cand.models[[1]] = BM
Cand.models[[2]] = OU
Cand.models[[3]] = Lambda
Cand.models[[4]] = EB
Cand.models[[5]] = OLS
Modnames = paste(c("BM", "OU", "Lambda", "EB", "OLS"), sep = " ")
aictab(cand.set = Cand.models, modnames = Modnames, sort = T)
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## Lambda 4 73.83 0.00 0.60 0.60 -32.73
## OU 3 75.29 1.46 0.29 0.89 -34.54
## OLS 3 77.15 3.32 0.11 1.00 -35.46
## EB 3 100.99 27.16 0.00 1.00 -47.39
## BM 3 155.40 81.57 0.00 1.00 -74.59
summary(Lambda)
## Generalized least squares fit by maximum likelihood
## Model: log(Eggshell_thickness) ~ log(Egg_mass)
## Data: datahard
## AIC BIC logLik
## 73.45572 84.32971 -32.72786
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.3282094
## Variance function:
## Structure: fixed weights
## Formula: ~Wh
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 4.321066 0.12577804 34.35469 0
## log(Egg_mass) 0.439690 0.01721966 25.53422 0
##
## Correlation:
## (Intr)
## log(Egg_mass) -0.483
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.87276697 -0.89716270 -0.55727271 -0.03185514 2.23701950
##
## Residual standard error: 0.02058541
## Degrees of freedom: 112 total; 110 residual
dataS<-read.table("Dataset2-amniotes.txt", header=T)
datahard<-dataS[-c(50,86:88,93:95,115,118,120:145,147,149),]
datacomp<-comparative.data(phy=treehard, data=datahard, names.col="Taxon")
pgls<-pgls(log(Eggshell_thickness)~log(Egg_mass), data=datacomp, lambda="ML")
summary(pgls)
##
## Call:
## pgls(formula = log(Eggshell_thickness) ~ log(Egg_mass), data = datacomp,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.06443 -0.01539 -0.00283 0.01159 0.05560
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.328
## lower bound : 0.000, p = 0.0012278
## upper bound : 1.000, p = < 2.22e-16
## 95.0% CI : (0.098, 0.597)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.334504 0.128404 33.757 < 2.2e-16 ***
## log(Egg_mass) 0.436835 0.017409 25.092 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02103 on 110 degrees of freedom
## Multiple R-squared: 0.8513, Adjusted R-squared: 0.8499
## F-statistic: 629.6 on 1 and 110 DF, p-value: < 2.2e-16
datasoft<-dataS[c(50,86:88,93,94,115,118,120:145,147,149),]
treesoft<-drop.tip(treeS, setdiff(treeS$tip.label, datasoft$Taxon))
dataS<-read.table("Dataset2-amniotes.txt", header=T, row.names="Taxon")
datasoft<-dataS[c(50,86:88,93,94,115,118,120:145,147,149),]
Ws<-diag(vcv.phylo(treesoft))
alpha <- seq(0, 1, 0.1)
fit <- list()
form <- log(Eggshell_thickness)~log(Egg_mass)
for (i in seq_along(alpha)) {
cor <- corMartins(alpha[i], phy = treesoft, fixed = T)
fit[[i]] <- gls(form, correlation = cor, data = datasoft, na.action=na.exclude, weights=varFixed(~Ws), method = "ML")
}
sapply(fit, logLik)
## [1] Inf -31.89120 -31.76165 -31.77510 -31.78728 -31.79212 -31.79376
## [8] -31.79428 -31.79444 -31.79448 -31.79450
g <- seq(0.1, 1, 0.1)
fit <- list()
form <- log(Eggshell_thickness)~log(Egg_mass)
for (i in seq_along(g)) {
cor <- corBlomberg(g[i], phy = treesoft, fixed = T)
fit[[i]] <- gls(form, correlation = cor, data = datasoft, na.action=na.exclude, weights=varFixed(~Ws), method = "ML")
}
sapply(fit, logLik)
## [1] -31.47980 -33.92230 -36.02601 -37.59877 -38.79144 -39.72307 -40.47140
## [8] -41.08692 -41.60329 -42.04372
BM<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datasoft, correlation=corBrownian(phy=treesoft), weights=varFixed(~Ws), method="ML")
OU<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datasoft, correlation=corMartins(1, phy=treesoft), weights=varFixed(~Ws), method="ML")
Lambda<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datasoft, correlation=corPagel(0.2, phy=treesoft), weights=varFixed(~Ws), method="ML")
EB<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datasoft, correlation=corBlomberg(0.1, phy=treesoft, fixed=T), weights=varFixed(~Ws), method="ML")
OLS<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datasoft, method="ML")
Cand.models = list()
Cand.models[[1]] = BM
Cand.models[[2]] = OU
Cand.models[[3]] = Lambda
Cand.models[[4]] = EB
Cand.models[[5]] = OLS
Modnames = paste(c("BM", "OU", "Lambda", "EB", "OLS"), sep = " ")
aictab(cand.set = Cand.models, modnames = Modnames, sort = T)
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## EB 3 69.71 0.00 0.51 0.51 -31.48
## Lambda 4 71.19 1.48 0.24 0.75 -30.95
## OLS 3 72.28 2.57 0.14 0.90 -32.77
## OU 4 72.88 3.17 0.10 1.00 -31.79
## BM 3 90.84 21.13 0.00 1.00 -42.04
summary(Lambda)
## Generalized least squares fit by maximum likelihood
## Model: log(Eggshell_thickness) ~ log(Egg_mass)
## Data: datasoft
## AIC BIC logLik
## 69.89509 76.22916 -30.94754
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.2416992
## Variance function:
## Structure: fixed weights
## Formula: ~Ws
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 2.1224777 0.21594179 9.828935 0
## log(Egg_mass) 0.3930565 0.05253711 7.481503 0
##
## Correlation:
## (Intr)
## log(Egg_mass) -0.691
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.2225960 -0.7969044 -0.3082419 0.5311175 1.5319272
##
## Residual standard error: 0.03652805
## Degrees of freedom: 36 total; 34 residual
dataS<-read.table("Dataset2-amniotes.txt", header=T)
datasoft<-dataS[c(50,86:88,93,94,115,118,120:145,147,149),]
datacomp<-comparative.data(phy=multi2di(treesoft), data=datasoft, names.col="Taxon")
pgls<-pgls(log(Eggshell_thickness)~log(Egg_mass), data=datacomp, lambda="ML")
summary(pgls)
##
## Call:
## pgls(formula = log(Eggshell_thickness) ~ log(Egg_mass), data = datacomp,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.103790 -0.018297 -0.004759 0.035618 0.068005
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.235
## lower bound : 0.000, p = 0.21201
## upper bound : 1.000, p = 2.3519e-06
## 95.0% CI : (NA, 0.751)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.113363 0.218395 9.6768 2.687e-11 ***
## log(Egg_mass) 0.390266 0.051037 7.6467 6.899e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03763 on 34 degrees of freedom
## Multiple R-squared: 0.6323, Adjusted R-squared: 0.6215
## F-statistic: 58.47 on 1 and 34 DF, p-value: 6.899e-09
dataS<-read.table("Dataset2-amniotes.txt", header=T)
ggplot(dataS, aes(log(Egg_mass), log(Eggshell_thickness), colour=Clade)) +
geom_point(size=4) +
xlab("ln egg mass (g)") +
ylab("ln calcareous layer thickness (µm)") +
geom_abline(intercept=4.321066, slope=0.439690, colour="turquoise", size=1.3) +
geom_abline(intercept=2.1224777, slope=0.3930565, colour="red", size=1.3) +
theme(panel.background = element_rect(fill="black")) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_colour_brewer("Group", palette="Set1")
dataS<-read.table("Dataset2-amniotes.txt", header=T); dataS<-dataS[-95,]
treeplot<-drop.tip(treeS, setdiff(treeS$tip.label, dataS$Taxon))
dataS<-read.table("Dataset2-amniotes.txt", header=T, row.names="Taxon"); dataS<-dataS[-95,]
dataplot=as.matrix(dataS[,3]); names(dataplot)<-rownames(dataS); colnames(dataplot)<-"CL_thickness"
dataplot<-log10(dataplot)
fit<-fastAnc(treeplot,dataplot,vars=T,CI=T) # Ancestral states for each node
obj<-contMap(treeplot,dataplot,plot=F)
plot(setMap(obj, colors=rev(brewer.pal(10, "Spectral"))), fsize=0.4,lwd=4)
plot(setMap(obj, colors=rev(brewer.pal(10, "Spectral"))), type="fan", fsize=0.4,lwd=4)
treeplot<-force.ultrametric(treeplot)
fit<-fastAnc(treeplot,dataplot,vars=T,CI=T) # Ancestral states for each node
obj<-contMap(treeplot,dataplot,plot=F)
plot(setMap(obj, colors=rev(brewer.pal(10, "Spectral"))), fsize=0.4,lwd=4)
dataS<-read.table("Dataset2-amniotes.txt", header=T)
dataS<-dataS[-95,]
treeS<-read.nexus("Amniotetree.nex")
treeS<-drop.tip(treeS, setdiff(treeS$tip.label, dataS$Taxon))
prismatic<-dataS[,5]; names(prismatic)<-dataS$Taxon
cols<-setNames(c("royalblue","red3"),levels(prismatic))
mtree<-make.simmap(treeS, prismatic, model="ER")
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## Absent Present
## Absent -0.0003949281 0.0003949281
## Present 0.0003949281 -0.0003949281
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Absent Present
## 0.5 0.5
## Done.
plot(mtree,cols,type="fan",fsize=0.5,ftype="i")
add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
y=0.8*par()$usr[3],fsize=0.8)
# with 1000 simulations of stochastic character maps from the data:
mtrees<-make.simmap(treeS, prismatic, model="ER", nsim=1000)
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## Absent Present
## Absent -0.0003949281 0.0003949281
## Present 0.0003949281 -0.0003949281
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Absent Present
## 0.5 0.5
## Done.
par(mfrow=c(10,10))
null<-sapply(mtrees,plotSimmap,colors=cols,lwd=1,ftype="off")
pd<-summary(mtrees); pd
## 1000 trees with a mapped discrete character with states:
## Absent, Present
##
## trees have 3.749 changes between states on average
##
## changes are of the following types:
## Absent,Present Present,Absent
## x->y 1.013 2.736
##
## mean total time spent in each state is:
## Absent Present total
## raw 2714.5230084 6746.7366863 9461.26
## prop 0.2869093 0.7130907 1.00
par(mfrow=c(1,1))
plot(pd,fsize=0.6,ftype="i",colors=cols,ylim=c(-2,Ntip(treeS)))
add.simmap.legend(colors=cols[2:1],prompt=FALSE,x=0,y=-2,vertical=FALSE)
# density map of the above
objprism<-densityMap(mtrees,states=levels(prismatic)[2:1],plot=FALSE)
## sorry - this might take a while; please be patient
n<-length(objprism$cols)
objprism$cols[1:n]<-colorRampPalette(c("royalblue","red3"))(n)
plot(objprism,size=c(0.6,1),fsize=0.6,lwd=4)
dataS<-read.table("Dataset2-amniotes.txt", header=T)
RTM<-dataS$Eggshell_thickness/dataS$Egg_mass
dataS<-cbind(dataS, RTM)
dataSboxplot<-dataS %>% filter(!Clade=="Test")
ggplot(dataSboxplot, aes(x=Clade, y=log(RTM), fill=Clade)) +
geom_boxplot() +
theme(panel.background = element_rect(fill="white")) +
theme(axis.text=element_text(size=10, angle = 90),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_fill_brewer(palette="Set1")
p<-ggplot(dataSboxplot, aes(x=Clade, y=log(RTM), fill=Clade)) +
geom_boxplot()
ggplot_build(p)$data
## [[1]]
## fill ymin lower middle upper ymax
## 1 #F8766D -4.1737717 -4.1737717 -4.1737717 -4.1737717 -4.1737717
## 2 #D39200 0.4296871 1.8265428 2.4120341 2.9394662 4.1207778
## 3 #93AA00 0.2636401 1.6216397 2.5162801 2.5794171 2.9038769
## 4 #00BA38 1.4714471 1.4714471 1.4714471 1.4714471 1.4714471
## 5 #00C19F 1.1282254 1.5449861 1.7324607 1.9034789 2.2137481
## 6 #00B9E3 4.4770380 4.4770380 4.4770380 4.4770380 4.4770380
## 7 #619CFF -2.1102132 -0.4793495 -0.1364957 0.9978001 1.9902104
## 8 #DB72FB 0.2876821 0.3729137 0.4581454 0.5433770 0.6286087
## 9 #FF61C3 -1.4307461 0.2392498 1.7701329 3.4510661 5.4647026
## outliers notchupper notchlower x flipped_aes y PANEL
## 1 -4.1737717 -4.1737717 1 FALSE -4.173772 1
## 2 -0.9350988, -0.6852781 2.6607121 2.1633561 2 FALSE NA 1
## 3 3.0207095 2.0118506 3 FALSE NA 1
## 4 1.4714471 1.4714471 4 FALSE 1.471447 1
## 5 1.9327199 1.5322015 5 FALSE NA 1
## 6 4.4770380 4.4770380 6 FALSE 4.477038 1
## 7 0.2372268 -0.5102183 7 FALSE NA 1
## 8 0.6485919 0.2676989 8 FALSE NA 1
## 9 2.5933530 0.9469128 9 FALSE NA 1
## group ymin_final ymax_final xmin xmax xid newx new_width weight colour size
## 1 1 -4.1737717 -4.1737717 0.625 1.375 1 1 0.75 1 grey20 0.5
## 2 2 -0.9350988 4.1207778 1.625 2.375 2 2 0.75 1 grey20 0.5
## 3 3 0.2636401 2.9038769 2.625 3.375 3 3 0.75 1 grey20 0.5
## 4 4 1.4714471 1.4714471 3.625 4.375 4 4 0.75 1 grey20 0.5
## 5 5 1.1282254 2.2137481 4.625 5.375 5 5 0.75 1 grey20 0.5
## 6 6 4.4770380 4.4770380 5.625 6.375 6 6 0.75 1 grey20 0.5
## 7 7 -2.1102132 1.9902104 6.625 7.375 7 7 0.75 1 grey20 0.5
## 8 8 0.2876821 0.6286087 7.625 8.375 8 8 0.75 1 grey20 0.5
## 9 9 -1.4307461 5.4647026 8.625 9.375 9 9 0.75 1 grey20 0.5
## alpha shape linetype
## 1 NA 19 solid
## 2 NA 19 solid
## 3 NA 19 solid
## 4 NA 19 solid
## 5 NA 19 solid
## 6 NA 19 solid
## 7 NA 19 solid
## 8 NA 19 solid
## 9 NA 19 solid